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Abstract 

RNA viruses exist in large intra-host populations which display great genotypic and phe- 
notypic diversity. We analyze a model of viral competition between two viruses infecting a 
constantly replenished cell pool, in which we assume a trade-off between the virus' colonization 
skills (cell killing ability or virulence) and its local competition skills (replication performance 
within coinfected cells). We characterize the conditions that allow for viral spread by means 
of the basic reproductive number and show that a local coexistence equilibrium exists, which 
is asymptotically stable. At this equilibrium, the less virulent competitor has a reproductive 
advantage over the more virulent colonizer. The equilibria at which one virus outcompetes the 
other one are unstable, i.e., a second virus is always able to permanently invade. One gen- 
eralization of the model is to consider multiple viral strains, each one displaying a different 
virulence. However, to account for the large phenotypic diversity in viral populations, we con- 
sider a continuous spectrum of virulences and present a continuum limit of this multiple viral 
strains model that describes the time evolution of an initial continuous distribution of viru- 
lence. We provide a proof of the existence of solutions of the model's equations and present 
numerical approximations of solutions for different initial distributions. Our simulations suggest 
that initial continuous distributions of virulence evolve towards a stationary distribution that is 
extremely skewed in favor of competitors. Consequently, collective virulence attenuation takes 
place. This finding may contribute to understanding the phenomenon of virulence attenuation, 
which has been reported in previous experimental studies. 

Keywords: SIR models of viral infection, Competition-colonization dynamics, RNA virus, Evolu- 
tion of virulence, Attenuation of virulency. 

1 Introduction 

RNA viruses are fast evolving pathogens that can adapt to continuously changing environments. 
Due to their error-prone replication, large population size, and high turnover, RNA virus popu- 
lations exist as quasispecies (Holland et al. (13), Domingo and Holland (8)). The viral mutant 
spectrum consists of many genetic variants which give rise to diverse phenotypes. This phenotypic 
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diversity is reflected in different traits, including the ability to kill host cells (herein referred to as 
virulence) . 

The concept of virulence has been used in various areas of the life sciences with different mean- 
ings. For instance, in evolutionary biology the virulence of a pathogen is defined as the fitness 
costs to the host that are induced by the pathogen. In epidemiology the term usually means the 
pathogen-induced host mortality. In more clinical settings virulence often refers to the severity of 
disease symptoms induced by a pathogen. In this article, we consider intra-host virus dynamics 
and use the term virulence to denote the cell killing rate of a virus infecting tissue (Herrera et al. 
(fl2l )). Thus, we apply the epidemiological meaning of the concept of virulence to the intra-host 
viral microepidemics. This use of the concept of virulence is related to the macroscopic or inter-host 
significations just explained, because, in general, the cell killing rate of a virus affects the course of 
infection, its inter-host repercussions as well as the mortality of the host. 

The evolution of virulence has been vastly studied in experimental and theoretical approaches 
in a variety of host-pathogen interactions and under diverse conditions or assumptions (see, among 
many others, Bull (0), May and Nowak (21), Lenski and May 0), Bonhoeffer et al. 0), Lenski 



18), Lenski and Levin (|19l ). Taylor et al. (|28l ). Haraguchi and Sasaki ([ill ). O'Keefe and Antonovics 

Different mechanisms have been proposed for the frequently observed coexistence of several 
viral strains. Quasispecies theory suggests that viral variants are maintained in a mutation-selection 
equilibrium (Eigen et al. (0)). In ecology, coexistence often results from spatially structured habitats 
(Tilman (29)). A trade-off between the ability of each individual to colonize unoccupied territory 
and to compete with others for the same habitat patch can result in coexistence of two strategies: 
competition and colonization. Competitors have an advantage when competing locally for resources, 
whereas colonizers are more successful in reaching new resources. 

Competition-colonization dynamics in an RNA virus infection have been recently demonstrated 
in vitro by Ojosnegros et al. (|25l ) in an experimental and theoretical approach. In this system, 
highly virulent viral strains play the role of colonizers, because they kill cells faster and thus 
replicate faster which allows faster spread and colonization of new cells. Local competition arises 
when two or more different viruses infect the same cell and compete for intracellular resources. 
Competitors manage to produce more offspring in a cell coinfected together with a colonizer and, 
at the same time, extend the cell killing time characteristic of a colonizer. This phenomenon is 
termed viral interference. 

The competition-colonization revolutionar y dy namics of two different viral strains in cell cul- 
ture have been described in Ojosnegros et al. (|25l ) by a modification of the basic model of virus 
dynamics (Nowak and May (|23l ). Perelson and Nelson (27|)). The model predicts that the outcome 
of viral competition for the cell monolayer depends on the initial density of viruses per cell. Under 
low initial density, colonizers produce more total offspring, whereas under high-density conditions, 
coinfection is more likely to occur and hence competitors have a selective advantage. This pre- 
diction was confirmed experimentally. Thus, it is plausible to think that low virulent viral strains 
(reduced cell killing ability, lack of colonization skills), which possess a reproductive advantage 
within coinfected cells (competition skills), will be naturally selected in a cell culture infection 
with multiple viral strains. In the same line of thought, a stronger relationship between both 
types of skills is conceivable, namely, a trade-off between virulence (colonization) and reproductive 
advantage within coinfected cells (competition). 

In the present article, we make a first step towards transferring the results obtained in vitro 
(Ojosnegros et al. (0)) to the in vivo situation under the assumption of a competition-colonization 
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trade-off. We extend the basic model of competition-colonization dynamics of two different viral 
strains in cell culture by replacing the finite cell monolayer with a constantly replenished pool of 
uninfected cells. Furthermore, to account for the competition-colonization trade-off postulated, we 
model the intracellular fitness of the viruses during coinfection as being inversely proportional to 
their respective virulence. In this relationship we archetypically capture the trade-off assumption. 

We present a rigorous analytical study of this model and demonstrate that it allows for local 
asymptotically stable coexistence of competitors and colonizers. Moreover, the less virulent com- 
petitor is shown to have a reproductive advantage reflected by its higher abundance and a higher 
number of cells infected with it at equilibrium. This counterintuitive result is in contrast to tra- 
ditional modeling approaches of viral competition, in which the equilibrium state is dominated by 
the most virulent strains (May and Nowak (|2ll)). 

A straightforward generalization of this model is to consider more than two viral variants. If we 
assume that different viral strains can only be distinguished through their virulences, the question 
arises as to how a distribution of virulences at the onset of infection is modified in the course of 
the infection by the competition-colonization dynamics. In this sense, the model describes the 
time evolution of a discrete distribution of virulences. While simulation results for a finite number 
of viral strains will be presented elsewhere (Ojosnegros et al., in preparation), herein we account 
for the very high phenotypic diversity of viral populations by considering the continuum limit of 
this multiple viral strains model. In this continuum limit, we consider a continuous spectrum of 
virulence values and associate a different viral strain with each virulence value. Thus, the time 
evolution of a continuous distribution of virulence is modeled. For this model's equations, we 
provide a proof of the existence of solutions and present numerical approximations of solutions for 
different initial distributions. 

Our simulations suggest that initial continuous distributions of virulence evolve towards a sta- 
tionary distribution that is extremely skewed in favor of the competitors. Thus, the model pre- 
dicts attenuation of the virus population and it might explain previous observations of suppres- 
sion of high- fitness mutants in various viral systems (de la Torre and Holland (|30l ). Novella et al. 
(0), Turner and Chao Jill), Bull et al. 0)). 

This article is organized as follows. In Section [2] we formulate the basic model of two competing 
viral strains and illustrate our model assumptions. In the Results Section a detailed analytical 
study of its equilibrium behavior is presented. In Subsection 13.51 the multiple- viral-strains model 
is introduced and in Subsection 13.61 we derive its continuum limit. The continuous- virulence model 
is analyzed both analytically and numerically. We close in Section 0] with discussing some of the 
model assumptions and consequences for the evolution of virulence. 



2 Formulation of the two-viral-strains model 

Our modeling approach is based on the basic SIR model of virus dynamics (Nowak and May 
I24I ) .Perelson and Nelson (j27h ). which we extend in order to model two different viral populations 
infecting constantly replenished tissue. We model singly infected and superinfected (doubly infected 
or coinfected) cells. Thus, the time evolution of the concentration of two competing viral strains, 
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uninfected and infected cells is described by the following system of ordinary differential equations 

x = A — dx — ftx{y\ + v 2 ) 

yi = /3xvi - f3y±v 2 - a\y\ 

y 2 = (3xv 2 - f3y 2 vi - a 2 y 2 (1) 

Z/12 = P(yiV2 + J/2«l) - min(ai, a 2 )yi2 

vi = Ka\y\ + cK min(ai, a 2 )y\ 2 — uv\ 

v 2 = Ka 2 y 2 + (1 - c)Kmm(ai,a 2 )yi 2 - uv 2 

The variable x models the concentration of uninfected cells with an external constant supply of 
new cells at rate A, dying at a rate d and being infected with efficiency (3. The variable v\ resp. v 2 
describes the concentration of strain 1 resp. strain 2. The variable yi represents the concentration of 
cells infected solely with strain i. These cells die and release viral offspring at rate cij, the virulence 
of strain i. The variable y\ 2 models the concentration of cells infected with both viral strains. These 
cells die and release viral offspring at rate min(ai,a2) (more on this below). Free virus of type i is 
produced at rate ki = Kai, where K is the burst size, and inactivated at rate u. The parameter 
c denotes the proportion of strain 1 produced at the burst of coinfected cells y± 2 . The state of the 
system at time t is denoted S(t) = (x(t), yx(t), y 2 (t), yi 2 {t), vi(t), v 2 {t)) T . 
We make the following general assumptions about the parameters: 

• All parameters are positive. 

• The efficiency with which strain 1 resp. strain 2 infects uninfected cells or singly infected cells 
(yi or y 2 ) is equal and denoted by j3. 

• The death rates of strain 1 resp. strain 2 are equal and denoted by u. 

Furthermore, based on the experimental results presented in Ojosnegros et al. jis). we assume: 

• The burst size of singly and coinfected cells is equal and denoted by K. 

• The death rate of coinfected cells is equal to the death rate of cells singly infected with the 
least virulent virus, 012 = min(ai, a 2 ). In other words, this rate is imposed by the least 
virulent viral strain. 

• To account for the competition-colonization trade-off postulated, we model the intracellular 
fitness of the viruses during coinfection as being inversely proportional to their respective 
virulence by setting c := a^ 1 /(a^ 1 + a^ 1 ). In this relationship we archetypically capture 
the assumption that, viruses compensate the lack of colonization skills with intracellular 
competition abilities and vice versa. 

Without loss of generality, we can assume a\ < a 2 . Thus, strain 1 is the competitor and strain 
2 is the colonizer and the last three equations of (H|) become 

2/12 = Pyiv 2 + /3y 2 vi - aiy 12 
v\ = Ka\y\ + cKa\y\ 2 — uv\ 
i) 2 = Ka 2 y 2 + (1 - o)Ka\y\ 2 - uv 2 
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In Ojosnegros et al. (|29). the model (P) with A = 0, d = 0, a\ < a 2 and an unconstrained (i.e. 
independent of ax) parameter c > 1/2, was introduced to describe two competing viral strains in 
cell culture. These special in vitro conditions of a fixed and limited amount of target cells allowed 
for an analytical treatment of the system in the large initial virus load limit which is not possible 
in the more general case of system ([!]) considered here. 



3 Results 

3.1 Establishing infection 

To identify the conditions on the parameters of the model that imply spread of at least one of the 
viral strains, we analyze the stability of the (obvious) equilibrium point 

s<°> = (xW$\$\i®A \<$ l >) T = (AM0,0,0,0,0f 

at which the infection dies out. The Jacobian matrix J of system ([I]) evaluated at is 



J(S<°>) 
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(1 - c)Kai 





-u j 



The eigenvalues of this matrix are 

-d(ai + u) + Aj 



-d, 



-ax, 



2d 



-d(a,i + u) — Aj 
2d 



1,2, 



where Aj := \J d 2 (ai — u) 2 + ^Ka-fiXd for i = 1,2. All eigenvalues are real given that all parameters 
are assumed to be positive. This equilibrium becomes unstable as soon as at least one eigenvalue 
is positive. This happens if and only if — d{a\ + u) + Ai > or — d(a2 + u) + A2 > which is 



equivalent to d[a\ 



u 



+ 4:Kai/3X > d(ai + u) 2 or d(a 2 - u) 2 + 4Ka 2 (3X > d(a 2 + u) 2 . The latter 



expression is in turn equivalent to Kj3\ > du. In other words, it is sufficient for viral spread that 

K/3X 



du 



> 1 



(2) 



For generic parameter values (the fine tuning K(3\ = du cannot be expected), this condition is also 
necessary, because K(3\ < du implies that is asymptotically stable. 

If we consider initial conditions in which yi(Q) = 0, 2/12(0) = 0, Vi(0) = and Vj(0) 7^ 0, where 
i,j E {l,2},i ^ j, the model reduces to a simple SIR model of single viral infection (Bonhoeffer 
et al. (3), Nowak and May (24), see also Korobeinikov (0) for global results) and we recognize the 
magnitude Rq as the well known basic reproductive number of the infection system. Condition ([2]) 
can also be expressed as M := Kf3\ — du > 0. The magnitude M turns out to be algebraically very 
helpful for our further analysis of the model. 

Having determined a condition on the parameter values that characterizes the event of viral 
spread, we next ask whether under these circumstances, the system admits a steady state in which 
both viral strains can coexist. The opposite steady state scenario would be that one of the viral 
strains outcompetes the other. To this end, we examine further fixed points of the system and their 
stability. 
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3.2 Further fixed points 



Performing algebraic manipulations we found the following non-trivial fixed points of the system 

/ u/K \ 



5 (D 



v4 1} 



i 

p 



5 (2) 



(2) 
I/l 
(2) 

y 2 

(2) 
2/12 
(2) 

V4 2 7 



V 



1 





M/( ai K) 



M/u 


u/K 


M/(a 2 K) 



M/u 



S* 



and 



5" 



y-2_ 







( u \ 


y* 




a 2 uM/(a 1 (M + u{a x + a 2 ))) 


y! 
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aiuAf/(a 2 (M + u(ai + a 2 ))) 


y* 2 




M 2 /(ai(M + u(ai +a 2 ))) 


«i 




a2KM/(u(ai + a 2 )) 


W 




\ aiKM/(u(ai + a 2 )) ) 


/ 




pK\/d 



1 



\ 



V 



(a 2 M - «dC)/((c - l)d(ai + a 2 )) 
-(aiM - ud£)/{cd(ai + a 2 )) 
-(aia 2 M + ud(((ai + a 2 )c - a 2 ))/(c(c - l)dai(ai + a 2 )) 

(K 
-(K 



where ^" is a root of the polynomial duZ 2 + (—c/3Ka 2 \ — fiK\a\c + f$K\a\ + da 2 u — a\ud)Z 
—aia 2 ud + f3K\a\a 2 in the indeterminate Z. Since we are only interested in real, non-negative 
solutions, the fixed point S~ can be disregarded, because it has at least one negative component for 
whatever real value £ might take (including £ = 0, which implies = —a\M/{cfiKd(a\ + a 2 ))). 
By the assumed positivity of the parameters and by M > (we assume that viral spread is 
possible) all other fixed points are (component-wise) non-negative and thus biologically meaningful. 
Summarizing, we have two equilibria, and , in which one of the viral strains outcompetes 
the other, and one coexistence equilibrium S*. The following subsections are devoted to studying 
their properties. 
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3.3 The coexistence equilibrium 



The Jacobian matrix J(S*) of system (TTJ) evaluated at the equilibrium point S* equals 

f-/3XK/u -u/K -u/K \ 

MM aiM n n ,,IK aiuM 

u(ai+a 2 ) u{ai+a 2 ) U1 U U U ' ^ a 1 K{M+u{a 1 +a 2 )) 

aiM n a2(M+u(a 1 +a2)) n amM n, I K 

u(ai+a 2 ) U u(oi+o 2 ) a 2 K(M+u(a 1 +a2)) U ' IX 

g a\M a 2 M a\uM a 2 uM 

u(ai+a 2 ) u(ai+a 2 ) 1 a 2 K (M +u(ai+a 2 )) a 1 K(M+u(a 1 +a 2 )) 



Kai *ogy -u 

V Ka * (S) "« J 

The characteristic polynomial of this matrix is too long to be displayed. Thus, we performed its 
analysis using computer algebra and symbolic computation. Under the premise of viral spread, i.e. 
M > 0, we used the computer algebra system Maple to show that all coefficients are positive. 
Furthermore, we used Maple^^ to construct Routh's table (Barnett and Siljak (1)) and verified 
that all entries in its first column are positive (see supplementary material). By Routh's criterion 
(Barnett and Siljak HI)), an ro °ts of the characteristic polynomial must lie strictly to the left of the 
imaginary axes. As a consequence, the equilibrium point S* in which both viral strains can coexist 
is a local asymptotically stable fixed point. This holds for all positive parameter values, provided 
M > 0. 

The peculiarity of this coexistence equilibrium is that the viral load at equilibrium of the 
competitor, v* = a<iM j {fiu(a\ + 02)), is proportional to the relative virulence of the colonizer. 
Similarly, the viral load at equilibrium of the colonizer, v\ = aiM/(f3u(a% + 02)), is proportional 
to the relative virulence of the competitor. Thus, the two-viral-strains system ([1]) not only allows 
for coexistence of both viral strains, but it confers the less virulent competitor a reproductive 
advantage over the more virulent colonizer. This discrepancy is reflected in the higher concentration 
of competitors, v\jv\ = 0,2/0,1 > 1, and the higher concentration of cells infected with competitors, 
2/1/2/2 = ( a 2/oi) 2 > 1, at equilibrium. The property that makes this phenomenon possible is the 
inverse proportionality between virulence and intracellular fitness during coinfection, as explained 
in more detail in the Discussion. 



3.4 Single viral strain equilibria 

The existence of a local asymptotically stable coexistence equilibrium suggests that a viral strain 
can bear the presence of another one. However, to fully address the question as to whether system 
(PQ) always allows for a second viral strain to invade tissue already infected with a different strain, 
we examine the stability of the equilibria and , in which one of the viral strains extrudes 
the other. 

The Jacobian matrix of system ([!]) evaluated at the point is 



J(S<V) 



M/u - d 











-u/K 


-u/K \ 


M/u 


-Ol 








u/K 


-M/( ai K) 








M 








u/K 








M/u 


-Ol 





M/( ai K) 





Ka x 





Kaia,2/{ai + a 2 ) 


— u 











Ka 2 


Ka\/(a\ + a 2 ) 





-u ) 
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Assuming viral spread, (M > 0), we showed that the leading coefficient of the characteristic 
polynomial of J(S^) is positive, whereas the independent coefficient is negative (see supplementary 
material). By Routh's criterion, at least one root of the characteristic polynomial has positive real 
part. Therefore, the equilibrium point in which the competitor outcompetes the colonizer is not 
stable. Analogously, we found that the equilibrium point in which the colonizer outcompetes 
the competitor is also unstable. Both statements hold true for all positive parameter values, 
provided that M > 0. 

However, it is worth mentioning that any trajectory S(t), for which one viral strain, say of type 
i, and all cells infected or coinfected with it have disappeared at some point in time s, would stay 
confined in the corresponding hyperplane Hi := {(x,yi,y2,yi2,vi,V2) T \ y% = yu = Vi = 0} for 
all t > s. As mentioned in Subsection 2, within the corresponding hyperplanes, the equilibria 
and become asymptotically stable, provided Rq > 1. Whether a trajectory starting outside Hi 
flows into the hyperplane or not remains to be analytically studied. Our simulations do not show 
any evidence for this type of behavior. 



3.5 Multiple-viral-strains model 



A straightforward generalization of our model (Tjjj that accounts for the experimentally observed 
diversity of viral populations is to consider more than two competing viral strains. This general- 
ization raises the question of how many viruses can coinfect a cell simultaneously. For example, 
the multiplicity of HIV-infected spleen cells has been reported between 1 and 8 with mean 3.2 
(Jung et al. (1 14)). However, for the sake of mathematical simplicity, we consider here the case in 
which at most two viruses can coinfect a cell. We assume that we can distinguish each of the viral 
strains via their corresponding virulences a%. As above, we assume inverse proportionality between 
virulence and intracellular fitness during coinfection. In other words, the proportion of strain i 

produced at the burst of cells yij coinfected with V{ and Vj is given by q := a" 1 / ^a" 1 + o-J 1 ^ j f° r 
all i,j € {l,...,n}, i < j , where n € N is the total number of viral strains modeled. Thus, the 
equations for the generalized model read: 



x 



Vi 



A — dx — fix ' 



3=1 



fixvi - fiyi 



am, i = !,■■■, n 



, 3=1 , 



Vij = P{VlVj + yjv{) -min(ai,aj)yij, I, j = 1, n and I < j 

( \ 



Kcnyi + a i 1 K 



^2 —i —Wi(l,j)min(ai,aj)yi 

. l,j a l +a 3 



uvi, i = 1, n 



J 



where Wi(l,j) = 1 if I = i or j = i, and otherwise Wi(l,j) = 0. The model does not explicitly account 
for the order of infection. Nevertheless, to increase the symmetry of the model and to simplify the 
notation, we will consider the order of infection events and separately model the populations yij 
and yji, where the order of the indices indicates the order of infection with the viral strains vi and 
Vj. With this notation, in general, yij ^ yji, and the variable yu in the two- viral-strains model ([TJ) 
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refers to yu + 2/21- To be consistent, we have to make sure that for I ^ j, the magnitude yij + yji 
obeys the corresponding equation, that is 



Vlj + Vji = P{yivj + yjvi) - mm(ai,a j )(yi j + yji) 
To ensure this, the equation for y^ becomes 

Vlj = PyiVj - min(a i ,a i )y/ j , I, j = 1, ...,n s.t. I / j 
Summarizing, we obtain the following model 



Ui 



A — dx 



(3xvi - f3y 

fiyivj 

Kcnyi + aT 1 



mm(ai,aj)yij, 
( 



E 



a-iyi 



, n 



1,3 = 1, -,n s.t. / / j 



. 3=1 "I 



+ a. 



—Kmm(ai,aj)(yij +yji 



(3) 



UVi, 



l,...,n 



Note that the magnitude ^ Vj(t) is the total viral population at any given point in time t. 

3=1 

Since the number of equations in this model grows quadratically with the number n of viral 
strains, it becomes rather involved to analyze it. In (Ojosnegros et al., in preparation) the results of 
numerical simulation for numerically tractable values of n are presented. Here, in compliance with 
the quasispecies view of viral populations, we devise a new approach to studying the evolution of 
virulence and consider the continuum limit of the multistrain model ([3]). In this continuum limit, 
we consider a continuous spectrum of virulence values and identify viral strains with virulence 
values. We call the resulting continuum limit the continuous-virulence model. In this model the 
viral quasispecies is represented by a time-dependent continuous distribution of virulence. Unlike 
the discrete multiple-viral-strains model, the continuum approach allows us to study the virulence 
distribution of diverse RNA virus populations in a manner independent of the number of different 
strain types. 



3.6 Continuous- virulence model 

We identify viral strains with their virulence a and denote by v(a,t) the density of viruses of 
type a at time t. If we consider an interval [01,02] C (0, 1) of possible virulences, then the initial 
distribution of virulences is defined by a continuos density function v(-,0) : M — > M which vanishes 
outside the interval [a±, 02]- The continuum limit of model ([3|) is therefore the following initial value 
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problem (also known as Cauchy problem): 

x(t) = A - dx(t) - px(t) J ' v(£,t)d£ 



ai 

' 0,2 



= 0x(t)v(a,t) - Py(a,t) \ / -«y(a.t) 



\ai 



dz 

— (a,b,t) = Py(a,t)v(b,t) - mm(a,b)z(a,b,t) (4) 
-g^(a,t) = Kay(a,t) + aT 1 K [ j _ 1 — min(a, b)(z(a, b, t) + z(b, a, t))db J —uv(a,t) 



Vh 



x(0) = x , y(^0) = j/ (0. ^,M,0) = zo(^M), «(e,0)=«o(0 

For every i G K>o, the function : [01,02] x [01,02] — >• M describes the density of coinfected 

cells with respect to the two-dimensional (Lebesgue) measure on E 2 . The value z(a,b,t) is only 
meaningful for our modeling purposes outside the diagonal a = b. Note that the exception j ^ i in 
the sums of the original discrete model can be neglected here because the values of a real valued 
function on a set of measure zero do not modify the value of the integral. If we assume that a 
solution of this initial value problem exists, then we can derive the following expressions: Setting 
V(t) := f® 2 v(£,t)d£ we have, for the first equation, x(t) = A — dx — j3x(t)V(t) and thus 



:(*) =1x0 + 1 Xe w ^dA 



x (t) = I x + I Xe w ^>d^J e~ w ® (5) 

where W[9) := J Q e d + pV{r)dr. The second equation can be solved as 

y(a,t)= ^y (a) + [3 J x(r)v(a, T)e Ua ^dr j e~ Ua ^ 

where U a (9) := f Q a + f3V(r)dT. Substituting the expression for x(t) into the expression for y(a,t) 
gives 

/ * / T \ \ 

o~Ua(t) 



y(a,t) = ^y (a)+(3 J v{a,r) ^x + j \e w ®<%j e u ^~ w ^d^j e~ u « 

)(o) +P J v(a,T) ^x + j \e w ®d$j e^'^drj e~ Ua{ ^ (6) 
The third equation yields 

z(a, b,t)= I zo(a, b) + J y(a, n)v(b, r])e m[ ^ a ' b)v dr] j e~ min(a ' b) * (7) 
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Substituting the expression for y(a,t) allows us to express z(a,b,t) + z(b,a,t) in terms of v(a,r) 
v(b,r), and integrals involving them as 



z(a,b,t) + z(b,a,t) = (z (a,b) + z (b,a) + R(a,b,t)) e -min ( a > fe )* 



where 



t 



R(a,b,t) := (5 / ((v(b, rj)y (a) + v(a, rj)y (b) + Q(a, b, 77)) e --(^-^)) ^ 







Q(a,M) := /3 y (v(MMa,T)e (a - d)T + i>(a, ^(6, r)e (f,-d)T ) |^x + J \e w ®d£j dr 

Finally, we substitute the expressions obtained for y(a,t) and z(a,b,t) into the differential 
equation for v(a,t) obtaining 



dv 
~dt 



— j min(a, b)(z (a,b) + z (b,a) + R(a,b,t))e- min ^ b)t db — uv(a,t) (8) 



+ cr l K 



A solution of the system necessarily has to fulfill this integro-partial differential equation for 
the function v(a,t). On the other hand, a solution of the latter equation that is continuous on 
[01,02] x R and partially differentiable with respect to t, and satisfies v(£,Q) = vq(£), allows for 
constructing a solution of the system (|4|) by means of substitution of v(a, t) into the expressions ([5]), 
([6]), and ([?]). Based on this idea, we show in Appendix [A] the existence of solutions of the system 

3.6.1 Simulation results 

In order to explore the dynamics of the continuous-virulence model, we numerically solved the 
Cauchy problem (j3J), as described in more detail in Appendix [Bj using typical parameter values 
given in Table 1. Because y, z, and v represent concentration densities, the units of y and v are 
[concentration] / [virulence] and the unit of z is [concentration] /[virulence] 2 . Given that the unit of 
virulence is [Time] -1 , y, z, and v are measured in units of [Time] / [Volume] , [Time] 2 /[Volume], and 
[Time] / [Volume] , respectively. The variable x represents a concentration and its unit is therefore 
[Volume] -1 . 

Starting from non infected tissue, i.e. xq = X/d, yo(£) = and Zq($, yu) = 0, we studied 
three different initial unnormalized continuous distributions of virulences given by the densities 
u(o,0) = v (a), a G [0.01,0.5] : 

1. A uniform distribution defined by vo(a) = 1 if a G [0.01, 0.5] and t>o(a) = otherwise 

2. A Gaussian distribution with mean value [i = 1/4 and standard deviation a = 1/40 given by 

VQ ( a ) = e -800(a-l/4)* 
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Parameter 


Description 


Value 


Units 


A 


Natural growth rate of uninfected population 


10 b 


(ml * h)~ L 


d 


Natural death rate of uninfected population 


0.05 


h- 1 


P 


Rate of infection 


5 • 10- 8 


ml/h 


K 


Burst size 


150 


Dimensionless 


u 


Clearance rate of free virus 


0.15 


h- 1 



Table 1: Parameters of the model of evolution of virulence during infection 

3. A mixture of two Gaussian distributions with mean values = 1/6, ^ = 1/3, equal standard 
deviations u\ = 02 = 1/300 and unequal weights Ai = 3/10, A2 = 7/10 given by vo(a) = 

(3/7) e -4500(a-l/6) 2 + (• 7/ / 10 ) e -4500(a-l/3) 2 

The time evolution of a flat initial virulence density is shown in snapshots in Figure 1. After a 
very short absorption phase, the density takes an exponential shape in favor of the most virulent 
part of the interval which is greatly amplified during the first eleven hours of simulation time. After 
approximately 130 hours a recession is observed, which by t = 140 h has already decreased the pop- 
ulations' densities by one order of magnitude. After approximately 190 hours a qualitative change 
takes place and the distribution starts loosing its exponential shape to become a non-symmetric 
unimodal distribution with mode around virulence = 0.257 at t = 240 h. This distribution starts 
traveling to the left, becomes narrower and more symmetric. By t = 600 h the distribution has 
moved further to the left and is now almost symmetric. At this point, the dynamics become sig- 
nificantly slower and we start observing an amplification effect. At t = 2534 h we observe a very 
narrow Gaussian distribution approaching the left boundary of the interval [0.01,0.5]. Here the 
distribution starts changing its shape to become an exponentially shaped distribution, this time in 
favor of low-virulence competitors. The changes become very slow and, at least numerically, the 
system seems to be reaching a stationary distribution, which is exponential and highly in favor of 
the smallest virulence values. (See Appendix [O for a similarly detailed description of Figures 2 and 
3.) 

Despite the prominent differences between the initial distributions, Figures 2 and 3 reveal similar 
qualitative properties of the dynamics, namely, a biphasic behavior comprising an initial phase in 
which the more virulent parts of the density are amplified, followed by a second phase in which 
the less virulent regions predominate. All three trajectories become very similar once the density 
becomes unimodal, although the time scales are significantly different, and all three reach very 
similar stationary distributions. In the case of the uniform and the Gaussian distribution, the 
stationary distributions reached differ only by small quantitative discrepancies within the same 
order of magnitude. 

For each of the three initial virulence densities vo(a), Figure 4 shows the time evolution of 
the expected value of the virulence, E(a)(t) = J^L^bv(b,t)/ \\v(b,t)\\ db, which is obtained by nor- 
malizing at each point in time, t, as the system of integro-partial differential equations is not 
norm-preserving. In this graph, we can clearly identify the two different regimes of the biphasic 
behavior described above. 
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x 10 
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io t=20700h 



0.02 0.04 0.06 0.08 0.1 
Virulence 



Figure 1: Time evolution of a uniform initial distribution of virulences. Each panel shows the shape 
of the density function at the point in time displayed in its titel. 



4 Discussion 

We have analyzed the evolution of virulence of an RNA viral quasispecies in which the cell killing 
capacity (what we call the virulence aj) of viruses is inversely related to the intracellular viral 
fitness within coinfected cells. In the case of two viral strains, this competition-colonization trade- 
off allows for stable coexistence of competitors and colonizers and each virus type can be invaded 
by the other, whenever the conditions for viral spread are given. These conditions do not depend 
on the particular virulence values a%. However, the population levels at the coexistence equilibrium 
do depend on the particular virulence values aj, and, as we saw above, this dependency is in favor 
of the least virulent viral strain. 

Generalizing this two-viral-strains model to multiple viral strains is conceptually straightfor- 
ward, but the resulting system of differential equations is difficult to study analytically. Moreover, 
the lack of accurate experimental measurements of the actual number of (in terms of virulence) 
different strains contained in a phenotypically diverse viral population limits the applicability of 
this modeling approach. Furthermore, the quadratic dependency of the number of equations on 
the number of strains n (recall Subsection 3.5) constrains the dimension of the models that can 
be numerically analyzed (Kryazhimskiy et al. (|17fl). We circumvented these issues by considering 
a continuous spectrum of virulence values and postulating a model that would describe the time 
evolution of a continuous distribution of virulences under the same type of competition-colonization 
trade-off. This model of continuous virulence is naturally derived as the continuum limit of the 
multiple-viral-strains model, thus providing a better modeling framework for the very high pheno- 
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Figure 2: Time evolution of a Gaussian mixture initial distribution of virulences. 



typic diversity of viral populations. While the model exhibits a complicated mathematical structure 
as an integro-differential Cauchy problem (Cushing da)), we were able to provide a simple proof 
of the existence of solutions. Having clarified the issue of existence of solutions, we solved this 
Cauchy problem numerically using typical parameter values. The discretization step size in the 
numerical scheme (see Appendix iBl) . clearly limits the accuracy of the numerical approximations. 
However, this numerical limitation does not represent a loss of modeling power, whereas, as ex- 
plained above, computational issues do limit the number of scenarios that can be modeled with the 
discrete multiple-strains model. 

Our simulation results indicate that the intra-host evolution of virulence is characterized by two 
phases. During the first phase, colonizers become more frequent and the average virulence of the 
population increases. In the second phase, the abundance of competitors increases and the mean 
population virulence decreases. Eventually, the virulence distribution appears to reach a steady 
state in which competitors dominate strongly over colonizers. 

In the two- viral-strains model ([I]), two major assumptions are the competition-colonization 
trade-off, c = a^ 1 /(a^ 1 + a^ 1 ), and the cell killing rate imposed by competitors in coinfected cells, 
i.e. the factor min(ai,a2) in the last three equations of model (pQ). The second assumption turns 
out to be not crucial, which can be seen as follows. If we replace the minimum in model (JT]) by the 
maximum and assume a\ < a 2 as before, then we obtain 

yi2 = (3(yiv 2 + j/2«i) - a 2 yi2 
v\ = Kaiyi + cKa 2 y\ 2 — uv\ 
v 2 = Ka 2 y 2 + (1 - c)Ka 2 yi2 - uv 2 
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Figure 3: Time evolution of a Gaussian mixture initial distribution of virulences. 

for the last three equations, while all others remain unchanged. Comparing this model to the 
original one we realize that strain 1 and strain 2 have interchanged their roles in the sense that 
the equation for strain 1 now has mixed virulence terms whereas the one for strain 2 has become 
homogeneous. In other words, we might as well write the maximum model as 



X 


= A — dx — f3x(vi + v 2 ) 


in 


= j3xv 2 - (3y 2 vi - a 2 y 2 


in 


= (3xvi - (3yiv 2 - ami 


yi2 


= Pyiv 2 + /3y 2 vi - a 2 yi 2 




= Ka 2 y 2 + cKa 2 yi 2 - uv 2 


i'l 


= Kami + (1 - c)Ka 2 yi 2 - uvi 



where c := a 2 /(oj + a 2 ). This model is equivalent to flTJ) with ai and a 2 having switched their 
roles. The coexistence equilibrium in the maximum model is 



(x*\ 




( u \ 


y*2 




aiuM/(a 2 (M + u{a 2 + 01))) 


y\ 


1 


a 2 uM/{ai(M + u(a 2 + ax))) 


* 

2/12 




M 2 /(a 2 (M + u(a 2 + ai))) 


v* 2 




aiMK/(u(a 2 + ax)) 






\ a 2 MK/(u(a 2 + ax)) ) 



at which competitors are more abundant than colonizers both as free virus, v\jv\ = 02/01 > 1, 



15 



Uniform initial distribution 




Gaussian initial distribution 

0.35 1 1 1 1 1 1 r- 



0.3 ■ 



Gaussian mixture initial distribution 



0.2 ■ 



0.1 ■ 



0.01 ■ 



0.01 ■ 



10 10 10 10 10 10 10 10 
time in hours (log scale) 




10 10 10 10 10 10 10 
time in hours (log scale) 




0.05 ■ 



0.01 ■ 



10 10 10 10 10 10 10 10 10 
time in hours (log scale) 



Figure 4: Time evolution of the expected value E(a) = J^^bvQ)^)/ \\v(b,t)\\ db of virulence for 
three different initial distributions. The crosses mark the value of E{a) at the time points (> 0) 
depicted in the corresponding evolution-of-distribution-figure (Figures 1-3). 



and inside cells, y\ly\ = ( a 2/°i) 2 > 1- The comparison of the minimum and the maximum models 
lets us conclude that the crucial property of both models that allows for a coexistence equilibrium 
in favor of competitors is not the cell killing rate imposed, but rather the inverse proportionality 
between virulence and intracellular fitness during coinfection. 

To investigate the consistency of our modeling approach, we compared the qualitative proper- 
ties of the dynamics displayed by the continuous- virulence model ([3J and the qualitative properties 
of the dynamics displayed be the discrete multiple-viral-strains model (J3j) (simulation results pre- 
sented in Ojosnegros et al., in preparation). Assuming that viral spread is possible (i.e., M > 0), 
we found that the basic biphasic feature, namely, an initial phase in which the more virulent strains 
colonize and expand, followed by a second phase in which the less virulent strains predominate, 
is present in the dynamics of both models. Also the shapes of the continuous distributions in 
the course of the simulations show similarity to the ones observed in the discrete case, as long as 
the discrete virulence values considered in the discrete model are uniformly distributed over the 
interval [01,02] C (0, 1) of possible virulences. These facts are not surprising, given that the sys- 
tem of ODEs solved to numerically approximate the solution of the continuous virulence system 
is structurally very similar to the equations of the discrete model, the only difference being the 
weights that appear in the Newton-Cotes formulas used to approximate the integrals v(£,t)d£ 
and f^ 2 a _ij|_ b _i min(a, b)(z(a, b, t) + z(b,a,t))db (see Appendix |B|) . Nevertheless, we found it very 
interesting that if we simulate the continuous- virulence model using parameters such that the con- 
dition for viral spread is not fulfilled (i.e., M < 0), the system evolves towards the zero density 
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function (results not explicitly shown). This outcome seems to be independent of the initial dis- 
tributions of virulence used. This result suggests that the condition for viral spread, which we 
originally derived for the two-viral-strains model, appears to be still correct in the continuum limit. 

On the other hand, the dynamics of the two models do show an important difference: Under 
viral spread conditions (i.e., M > 0), the stationary distributions reached in the continuous and the 
discrete case differ in that the stationary continuous distribution is extremely more positively skewed 
and only the very least virulent strains are represented. This finding has unexpected implications 
for the discrete model ([3]), which will be discussed elsewhere (Ojosnegros et al., in preparation). 

In conclusion, the two models studied in this article make two major predictions about the 
evolution of virulence under a competition-colonization trade-off. First, two viral strains with 
distinct virulence can coexist, and second, a viral population displaying a range of virulence values 
will be attenuated and evolve towards a population of many competitors and very few colonizers. 
Our model predictions differ from most prognoses based on previous models of the evolution of 
virulence, which often conclude that selection maximizes the basic reproductive number of the 
pathogen. This discrepancy is due to the following: 

1. One key feature of our approach is that coinfections are explicitly modeled. 

2. We do not assume that the viral strain with the highest individual cell killing performance 
dominates the events during coinfections. 

3. We assume a competition-colonization trade-off. 

These three assumptions combined have not been considered in previous modeling approaches. 
Under these premises, selection appears to favor low-virulence competitors, as long as uninfected 
cells are constantly replenished, but not unlimited. The attenuation property of the continuous- 
virulence model may also explain experimental observations of suppression of high-fitness viral 
mutants (colonizers), which might have been displaced by competitors (de la Torre and Holland 
(0), Novella et al. 0), Turner and Chao (0), Bull et al. 0)). 



In Ojosnegros et al. (|25l ) two foot-and-mouth disease viral strains were reported that had been 
isolated from a population undergoing viral passaging experiments. Measurement of cell killing 
rates, intracellular fitness, and other parameters suggested a competition-colonization trade-off. 
These experimental findings motivate the hypothesis that viruses can specialize either to improve 
their colonization skills by fast cell killing, or to improve their competitive intracellular reproductive 
success. In our modeling approach, we have implemented the competition-colonization trade-off 
using the algebraically simple relationship c = a^ 1 /(a^ 1 + a^ 1 ) which renders the mathematical 
analysis convenient. The actual dependency between virulence and intracellular fitness c is 
likely to be more complicated and to depend on additional parameters. It would be of biological 
interest to identify and to characterize virus populations with a competition-colonization trade-off 
and to establish the nature of the trade-off experimentally. In this article, our principal aim was to 
provide a rigorous mathematical analysis of a viral competition model incorporating a simple but 
archetypical instance of a competition-colonization trade-off. 
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Appendix 



«1 



A Existence of solutions of the continuous- virulence model's equa- 
tions 

In order to show that a solution of © exists, we consider solutions during a very short time 
span [ti,t2] C [0,oo) within which the values of V(t) and z(a,b,t) do not significantly change, i.e. 
V(t) « Vti '■= V(t\) and z(a, b, t) « 2^ (a, b) := 0(0, 6, ii) V t € [ii, £2]' Given that (jl]) is autonomous, 
we may as well consider the time interval [t\ = 0, £2]- Thus, W(9) « Jq d + pVodr = 6(d + /3Vo), 

t/o(9) « f°a + pVo<fr = 9(a + pV ) V # G [0,t 2 ] and / a _ 1 ^_ 1 mm{a,b)(z(a,b,t) + z(b,a,t))db 

Ol 

122 

J* a _ij|_ b _i min(a, 6)(^o(«, &) + £o(&j a))cZ6 =: 5*o( a ) V £ € [0, t 2 ]- With this, equation ([8]) becomes 

|(<M) = Ka ^o(a) + ^/ /3t,(a,r) [x Q + ^™ v ~ ^ e^d^j e ~**™ 

+a~ 1 KSo{a) — uv(a,t) 
where yo(a) = y(a,0), xq = x(0). Some algebra yields 

+a _1 K5o(a) — uv(a,t) 

t 

= -uv(a,t)+Kf3Xae- t( - a+ ^ [ (^\— e {a_d)T + — ^— e r ( a+/3Vb A v(a, r)dr 

J \d^ + d(3Vo d + pVo ) 



+Kay (a) + a~ 1 iC5'o(a) 

If we write e~ t ^ a ^ v °' as 1 + (—a — /3Vo)t + 0(i 2 ) and neglect terms of quadratic order we obtain 
for t sufficiently small 



dv 



[a,t) = -uv(a,t)+Kp\a(l + (-a-pV )t) J |^ e (-^ + _^ e ^Wo A u(a> r)d7 



+ifay (a) + a 1 KS , (a) 



Let us assume for a moment that a three times differentiable solution exists. Differentiation on 
both sides yields 

- -4(°-')-« + ^)/(^ eMT+ ^k^ + '' , '° ) )" (< '' T)dT 
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Again, differentiation on both sides gives 



£(M) ^ -^(M)-^A«(« + «( 5 4|^e(.-* + ^L_e* + «))„(<., t ) 
+* PHI + (-a - W , (3^^ + S^W*"'") * ( * " 

W ! + <-. - f,V M (fcf^» + |±f£e«« + «>) „(. t ) 
Summarizing we have 

f£ (a , t) + „^ (o , t) _ ^ + ( _ a _ mt) (^|^ <->< + j^-/^) £<• 
= - 2 ^° + w {-^rkv/^ + s+W*™) »(«■ ') 

For each a S [01,02] the latter equation is a third order ordinary differential equation with 
variable coefficients. Using standard Lipschitz-continuity arguments (see, for instance, Section 4.3 
in Konigsberger (fl~5h ) it can be shown that for each o € [01, 02] the initial value problem together 
with v(a,0) = vo(a), |j| (a, 0) = -muo(o) + Kayo(a) + a~ 1 KS'o(a) and 

f^r (a, 0) = -u (-wuo(a) + Kay (a) + a _1 KS (a)) + K(3Xa ( ^^v Q + H+|vb ) u o(a) (which we as- 
sume to be continuous functions of a) must have a unique solution defined on some interval 
[0, T{\ C [0,oo) of positive length T\ G R+. Given that the coefficients of Q are continuous 
functions of a (which can be interpreted as a parameter in the ODE © in the framework of a sen- 
sitivity analysis) all the solutions v(a,t) must be continuous on [01,02] x [0, T\] (see, for instance, 
Theorem 6.1 in Epperson (fiol ) and also Subsection 3.1.1 in Deuflhard and Bornemann (0))- This 
family of solutions allows us to construct a local solution (defined on [01,02] x [0, T\\) of (|4|) by 
means of substitution of v(a, t) into the expressions ([5]), ([6]) and ([7|). The procedure can be now 
repeated for a short time interval starting at T\ yielding the next local solution. A global solution 
can be obtained by patching together the local solutions and letting the Tj — > 0. 

B Numerical solution of the continuous- virulence model's equa- 
tions 

To solve the system of equations dH numerically, we discretized the "virulence-space" with an 
equidistant grid G n ([0.01, 0.5]) and approximated the integrals v(£, t)d£ ~ J2jeG n ([o 01 5]) lj v {ji 
and lai a-tlb-t ™in(a,b)(z(a,b,t) + z(b,a,t))db « EjeG„([o.oi,o.5])7j a -i+j-i min(a, j)(z(a, j, t) + 
z(j,a,t)) using a Newton-Cotes quadrature formula of seventh order (with weights jj ). After this 
discretization step, we obtain for each pair (a,b) € (G n ([0.01, 0.5]) x G n ([0.01, 0.5])) the following 
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system of ordinary differential equations 



x(t) 



\-dx{t)-Px(t) Yl TjwO'.t) 

jeG„([0.01,0.5]) 




f3x(t)v(a, t) - (3y(a, t) ^ jjv(j, t) - ay{a, t) 
jeGn ([0.01,0.5]) 




a 



1 



-i 



min(a, j)(z(a,j,t) +z(j,a,t)) - uv(, 



The resulting system of coupled ordinary differential equations is solved numerically. 

C Detailed description of Figures 2 and 3 

The time evolution of an initial Gaussian distribution of virulences is shown in snapshots in Figure 
2. Initially we observe amplification, followed by even stronger amplification in combination with 
slight right-shift. At t = 20.2 h this phase reaches a peak. A minor recession follows and left- 
shifting sets on. At t = 700 h left shifting has progressed and moderate amplification sets on again. 
As left-shifting continues the distribution becomes narrower and narrower. The dynamics become 
slower as the left traveling peak approaches the left boundary of the interval [0.01,0.5]. Once it 
reaches the boundary, the distribution starts loosing its Gaussian shape to become an exponentially 
shaped distribution in favor of the least virulent part. The changes become very slow and, at least 
numerically, the system seems to be reaching a stationary distribution, which is exponential and 
highly in favor of the smallest virulence values. 

The time evolution of an initial Gaussian mixture distribution of virulences is shown in snapshots 
in Figure 3. Initially we observe amplification, with clear advantage for the most virulent component 
of the Gaussian mixture. At t = 17.19 h this component reaches a peak, while the least virulent 
component continues growing (order of magnitude 10 7 , thus not visible on the scale depicted). 
After this point, recession sets on. At t = 22.23 h both components are shrinking, the least virulent 
one only slightly (not visible). After approximately 14 hours the least virulent component starts 
growing again, while the most virulent one continues shrinking. By t = 178.15 h we can already 
see the least virulent component reach the order of magnitude of 10 7 . By t = 446.52 the most 
virulent component is about to fall below this order of magnitude. The least virulent component 
continues growing while the most virulent one declines further. At t = 9600 h the most virulent 
component has vanished (numerical order of magnitude 10~ 104 ) and the least virulent one starts 
shifting to the left while further amplified. At t = 37600 h it becomes visible that the least virulent 
component starts also narrowing, while growing and left-shifting continues. The dynamics remain 
qualitatively the same, become significantly slower though. The changes become very slow and, 
at least numerically, the system seems to be reaching a stationary distribution with mode 0.0393. 
This value is approximately 38 x a\ left from /ii. 
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